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This paper provides a detailed introductory description of Subset Simulation, an advanced stochas¬ 
tic simulation method for estimation of small probabilities of rare failure events. A simple and 
intuitive derivation of the method is given along with the discussion on its implementation. The 
method is illustrated with several easy-to-understand examples. For demonstration purposes, the 
MATLAB code for the considered examples is provided. The reader is assumed to be familiar only 
with elementary probability theory and statistics. 
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I. INTRODUCTION 

Subset Simulation (SS) is an efficient and elegant 
method for simulating rare events and estimating the 
corresponding small tail probabilities. The method was 
originally developed by Siu-Kui Au and James Beck in 
the already classical paper [1] for estimation of structural 
reliability of complex civil engineering systems such as 
tall buildings and bridges at risk from earthquakes. The 
method turned out to be so powerful and general that 
over the last decade, SS has been successfully applied to 
reliability problems in geotechnical, aerospace, fire, and 
nuclear engineering. Moreover, the idea of SS proved 
to be useful not only in reliability analysis but also in 
other problems associated with general engineering sys¬ 
tems, such as sensitivity analysis, design optimization, 
and uncertainty quantification. As of October 2013, ac¬ 
cording to the Web of Knowledge database, the original 
SS paper [T] received more than 250 citations that indi¬ 
cates the high impact of the Subset Simulation method 
on the engineering research community. 

Subset Simulation is essentially based on two different 
ideas: conceptual and technical. The conceptual idea is 
to decompose the rare event F into a sequence of pro¬ 
gressively “less-rare” nested events, 

F = C En-i C ... C Fi, (1) 

where Fi is a relatively frequent event. For example, 
suppose that F represents the event of getting exactly m 
heads when flipping a fair coin m times. If m is large, 
then F is a rare event. To decompose F into a sequence 
0 . let us define Ffc to be the event of getting exactly k 
heads in the first k flips, where k = 1,... ,m. The smaller 
k, the less rare the corresponding event F^; and Fi — 
getting heads in the first flip — is relatively frequent. 

Given a sequence of subsets 0 - the small probability 
P(F) of the rare event F can then be represented as a 
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product of larger probabilities as follows: 

P(F) = P(F™) 

= P(F^) _ 

^ '^P(Fi) P(F2) ■ ■ ■ P(F^_2) P(F„-i) ^ ^ 

= P(Fi) • P(F2 |Fi) • ... • P(F„|F„_i), 

where P(Ffc|Ffc_i) = P(Ffc)/P(F/j_i) denotes the condi¬ 
tional probability of event F^ given the occurrence of 
event F^-i, for k = 2,...,m. In the coin example, 
P(Fi) = 1/2, all conditional probabilities P(Ffe|F/j_i) = 
1/2, and the probability of the rare event P(F) = 1/2"*. 

Unlike the coin example, in real applications, it is of¬ 
ten not obvious how to decompose the rare event into a 
sequence p]) and how to compute all conditional proba¬ 
bilities in^^. In Subset Simulation, the “sequencing” 
of the rare event is done adaptively as the algorithm 
proceeds. This is achieved by employing Markov chain 
Monte Carlo, an advanced simulation technique, which 
constitutes the second ~ technical - idea behind SS. Fi¬ 
nally, all conditional probabilities are automatically ob¬ 
tained as a by-product of the adaptive sequencing. 

The main goals of this paper are: (a) to provide a de¬ 
tailed exposition of Subset Simulation at an introductory 
level; (b) to give a simple derivation of the method and 
discuss its implementation; and (c) to illustrate SS with 
intuitive examples. Although the scope of SS is much 
wider, in this paper the method is described in the con¬ 
text of engineering reliability estimation, the problem SS 
was originally developed for in 

The rest of the paper is organized as follows. Sec¬ 
tion describes the engineering reliability problem and 
explains why this problem is computationally challeng¬ 
ing. Section |ni] discusses how the direct Monte Carlo 
method can be used for engineering reliability estimation 
and why it is often inefficient. In Section [TV] a necessary 
preprocessing step which is often used by many reliability 
methods is briefly discussed. Section |V] is the core of the 
paper, where the SS method is explained. Illustrative ex¬ 
amples are considered in Section |VI| For demonstration 
purposes, the MATLAB code for the considered exam¬ 
ples is provided in Section |VII| Section |vnT] concludes 
the paper with a brief summary. 
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II. ENGINEERING RELIABILITY PROBLEM 

One of the most important and computationally chal¬ 
lenging problems in reliability engineering is to estimate 
the probability of failure for a system, that is, the proba¬ 
bility of unacceptable system performance. The behavior 
of the system can be described by a response variable y, 
which may represent, for example, the roof displacement 
or the largest interstory drift. The response variable de¬ 
pends on input variables x = {xi,...,Xd), also called 
basic variables, which may represent geometry, material 
properties, and loads, 

y = g{xi,...,Xd), (3) 

where g{x) is called the performance function. The per¬ 
formance of the system is measured by comparison of the 
response y with a specified critical value y*: ii y < y*, 
then the system is safe; ii y > y*, then the system has 
failed. This failure criterion allows to define the failure 
domain F in the input x-space as follows: 


Moreover, traditional numerical integration is also gen¬ 
erally not applicable. In this approach, the d-dimensional 
input x-space is partitioned into a union of disjoint hy¬ 
percubes, □ i,...,nAr. For each hypercube a “rep¬ 
resentative” point is chosen inside that hypercube, 
S Di- The integral in ^ is then approximated by 
the following sum: 

PF^ 7r(a;(*^)vol(ni), (6) 

where vol(ni) denotes the volume of and summation 
is taken over all failure points x^’''^. Since it is not known 
in advance whether a given point is a failure point or 
not (the failure domain F is not known explicitly), to 
compute the sum in <© , the failure criterion must 
be checked for all x^'^\ Therefore, the approximation ^ 
becomes 

N 

PF ^'^lF{x^'‘^)F{x^'‘'^)vol{n,), (7) 

2=1 


F = {x : g(x) > y*}. 


(4) where If(x) stands for the indicator function, i.e.. 


In other words, the failure domain is a set of values of 
input variables that lead to unacceptance system per¬ 
formance, namely, to the exceedance of some prescribed 
critical threshold y *, which may represent the maximum 
permissible roof displacement, maximum permissible in¬ 
terstory drift, etc. 

Engineering systems are complex systems, where com¬ 
plexity, in particular, means that the information about 
the system (its geometric and material properties) and its 
environment (loads) is never complete. Therefore, there 
are always uncertainties in the values of input variables 
X. To account for these uncertainties, the input vari¬ 
ables are modeled as random variables whose marginal 
distributions are usually obtained from test data, expert 
opinion, or from literature. Let 7r(x) denote the join 
probability density function (PDF) for x. The random¬ 
ness in the input variables is propagated through (|^ into 
the response variable y, which makes the failure event 
{a: € F} = {y > y*} also random. The engineering re¬ 
liability problem is then to compute the probability of 
failure pp, given by the following expression: 


If{x) 


1, if a; G F, 
0, if a; ^ F. 


( 8 ) 


If n denotes the number of intervals each dimension of 
the input space is partitioned into, then the total number 
of terms in Q is iV = n‘^. Therefore, the computational 
effort of numerical integration grows exponentially with 
the number of dimensions d. In engineering reliability 
problems, the dimension of the input space is typically 
very large (e.g., when the stochastic load time history is 
discretized in time). For example, d ~ 10^ is not un¬ 
usual in the reliability literature. This makes numerical 
integration computationally infeasible. 

Over the past few decades, many different methods for 
solving the engineering reliability problem ([^ have been 
developed. In general, the proposed reliability methods 
can be classified into three categories, namely: 

(a) Analytic methods are based on the Taylor-series ex¬ 
pansion of the performance function, e.g. the First- 
Order Reliability Method (FORM) and the Second- 
Order Reliability Method (SORM) [7l fT^ 1^ . 


Pf = P(a: GF)= Tr{x)dx. (5) 

Jf 

The behavior of complex systems, such as tall build¬ 
ings and bridges, is represented by a complex model (|^. 
In this context, complexity means that the performance 
function g{x), which defines the integration region F in 
([^, is not explicitly known. The evaluation of g{x) for 
any x is often time-consuming and usually done by the 
finite element method (FEM), one of the most important 
numerical tools for computation of the response of engi¬ 
neering systems. Thus, it is usually impossible to evalu¬ 
ate the integral in (§ analytically because the integration 
region, the failure domain F, is not known explicitly. 


(b) Surrogate methods are based on a functional sur¬ 
rogate of the performance function, e.g. the Re¬ 
sponse Surface Method (RSM) [UlTnillHl, Neural 
Networks [5S], Support Vector Machines [T^], and 
other methods [13]. 

(c) Monte Carlo simulation methods, among which 
are Importance Sampling [9|, Importance Sam¬ 
pling using Elementary Events [5], Radial-based 
Importance Sampling m, Adaptive Linked Im¬ 
portance Sampling [16], Directional Simulation |7], 
Line Sampling m, Auxiliary Domain Method m. 
Horseracing Simulation [3D], and Subset Simula¬ 
tion [T], 
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Subset Simulation is thus a reliability method which is 
based on (advanced) Monte Carlo simulation. 


III. THE DIRECT MONTE CARLO METHOD 

The Monte Carlo method, referred in this paper as Di¬ 
rect Monte Carlo (DMC), is a statistical sampling tech¬ 
nique that have been originally developed by Stan Ulam, 
John von Neumann, Nick Metropolis (who actually sug¬ 
gested the name “Monte Carlo” m), and their collab¬ 
orators for solving the problem of neutron diffusion and 
other problems in mathematical physics [22]. From a 
mathematical point of view, DMC allows to estimate the 
expected value of a quantity of interest. More specifi¬ 
cally, suppose the goal is to evaluate E^[/i(a;)], that is an 
expectation of a function h : X ^ M. with respect to the 
PDF 


E 7 r[/i(a;)] = / h{x)Tr{x)dx. (9) 

Jx 

The idea behind DMC is a straightforward application of 
the law of large numbers that states that if x^^\ ... 

are i.i.d. (independent and identically distributed) from 
the PDF 7r(a:), then the empirical average ^ h(a;^*^) 

converges to the true value E^[/i(x)] as N goes to -l-oo. 
Therefore, if the number of samples N is large enough, 
then E,r[h(a;)] can be accurately estimated by the corre¬ 
sponding empirical average: 


(10) 

The relevance of DMC to the reliability problem ([^ 
follows from a simple observation that the failure proba¬ 
bility pj’ can be written as an expectation of the indicator 
function (|^, namely. 


the fact that a;0) 


7r(x) and (111, 


E[p^“°] = E 


N 




1 ^ 
^ 2=1 


(13) 


The main advantage of DMC over numerical integra¬ 
tion is that its accuracy does not depend on the dimen¬ 
sion d of the input space. In reliability analysis, the 
standard measure of accuracy of an unbiased estimate 
Pf of the failure probability is its coefficient of vari¬ 
ation (c.o.v.) S(pf), which is defined as the ratio of 
the standa rd dev iation to the expected value of pp, i.e., 
Hpf) = \/V[pf]/E[pf], where V denotes the variance. 
The smaller the c.o.v. 6{pf), the more accurate the esti¬ 
mate pf is. It is straightforward to calculate the variance 
of the DMC estimate: 


V[p^“^] = V 




N 




1 

1 


2 = 1 
N 


^(E[/F(xW)2]-E[/F(a:«)]2 


2 = 1 
N 


= ]^J2(pp-Pf) = 


2 = 1 


Pf(I - Pf) 

N 


(14) 


Here, the identity If{x)^ = If{x) was used. Using (13) 
and (14), the c.o.v. of the DMC estimate can be calcu¬ 
lated: 


Pf = / Tr{x)dx = / lFix)Tr{x)dx = EjrllFix)], (11) 
Jf Jx 

where X denotes the entire input x-space. Therefore, 
the failure probability can be estimated using the DMC 
method (10) as follows: 


1 ^ 

PffspT- = -Y,If{x^% (12) 

2=1 

where are i.i.d. samples from 7r(a:). 

The DMC estimate of pf is thus just the ratio of the 
total number of failure samples If{x^'^'1), i.e., sam¬ 

ples that produce system failure according to the system 
model, to the total number of samples, N. Note that 
p™'^ is an unbiased random estimate of the failure prob¬ 
ability, that is, on average, p™^ equals to pp- Mathe¬ 
matically, this means that E[p™°] = pp- Indeed, using 


S{pT^ 


) 


E[p°“°] 


1 -PF 

NpF 


(15) 


This result shows that (5(p^“°) depends only on the fail¬ 
ure probability pp and the total number of samples N, 
and does not depend on the dimension d of the input 
space. Therefore, unlike numerical integration, the DMC 
method does not suffer from the “curse of dimensional¬ 
ity”, i.e. from an exponential increase in volume associ¬ 
ated with adding extra dimensions, and is able to handle 
problems of high dimension. 

Nevertheless, DMC has a serious drawback: it is inef¬ 
ficient in estimating small failure probabilities. For typ¬ 
ical engineering reliability problems, the failure proba¬ 
bility Pf is very small, pp <?; 1. In other words, the 
system is usually assumed to be designed properly, so 
that its failure is a rare event. In the reliability litera¬ 
ture, Pf 10“^ — I0“® have been considered. If pp is 
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very small, then it follows from (151 that 


where (/>(•) denotes the standard Gaussian PDF, 


s{prn 


Npp 


(16) 


This means that the number of samples N needed to 
achieve an acceptable level of accuracy is inverse propor¬ 
tional to pf, and therefore very large, N oc I/pf ^ 1- 
For example, \i pf = 10“^ and the c.o.v. of 10% is de¬ 
sirable, then N = 10® samples are required. Note, how¬ 
ever, that each evaluation of If^x^’''^), i = 1,...,1V, in 
(121 requires a system analysis to be performed to check 
whether the sample is a failure sample. As it has been 
already mentioned in Section |llj the computation effort 
for the system analysis, i.e., computation of the perfor¬ 
mance function g{x), is significant (usually involves the 
FEM method). As a result, the DMC method becomes 
excessively costly and practically inapplicable for reliabil¬ 
ity analysis. This deficiency of DMC has motivated re¬ 
search to develop more advanced simulation algorithms 
for efficient estimation of small failure probabilities in 
high dimensions. 


Remark 1. It is important to highlight, however, that 
even though DMC cannot be routinely used for reliability 
problems (too expensive), it is a very robust method, and 
it is often used as a check on other reliability methods. 


IV. PREPROCESSING: TRANSFORMATION 
OF INPUT VARIABLES 


Many reliability methods, including Subset Simula¬ 
tion, assume that the input variables x are independent. 
This assumption, however, is not a limitation, since in 
simulation one always starts from independent variables 
to generate dependent “physical” variables. Further¬ 
more, for convenience, it is often assumed that x are i.i.d. 
Gaussian. If this is not the case, a “preprocessing” step 
that transforms x to i.i.d. Gaussian variables z must be 
undertaken. The transformation form x to z can be per¬ 
formed in several ways depending on the available infor¬ 
mation about the input variables. In the simplest case, 
when X are independent Gaussians, Xk A/'(-|/rf;, cr^), 
where pk and are respectively the mean and variance 
of Xk, the necessary transformation is standardization: 


Zk 


^k 


(17) 


In other cases, more general techniques should be used, 
such as the Rosenblatt transformation m and the Nataf 
transformation [24j . To avoid introduction of additional 
notation, hereinafter, it is assumed without loss of gen¬ 
erality that the vector x has been already transformed 
and it follows the standard multivariate Gaussian distri¬ 
bution, 


(t>{x) 


1 




(19) 


V. THE SUBSET SIMULATION METHOD 

Unlike Direct Monte Carlo, where all computational 
resources are directly spent on sampling the input space, 
x^^\ ..., x^^'> ^ 7r(-), and computing the values of the 
performance function g{x ^^'>),... ,g{x^^'>), Subset Simu¬ 
lation first “probes” the input space X by generating a 
relatively small number of i.i.d samples Xq^\ ... ,Xq^'^ ^ 
7r(x), n < N, and computing the corresponding system 
responses = g(x^g^),...,y!f"^ = g(x[|"^). Here, the 
subscript 0 indicates the 0*^ stage of the algorithm. Since 
A is a rare event and n is relatively small, it is very likely 
that none of the samples Xq^\ ..., Xq”^ belongs to F, that 
is 2/o < 2 /* for alH = 1,..., n. Nevertheless, these Monte 
Carlo samples contain some useful information about the 
failure domain that can be utilized. To keep the notation 
simple, assume that yg^\ ..., ?/q"^ are arranged in the de¬ 
creasing order, i.e. > . ■. > (if i® always possible 
to achieve this by renumbering Xg^\ ..., Xq"^ if needed). 
Then, Xg^^ and Xg"^ are, respectively, the closest to fail¬ 
ure and the safest samples among Xg^^ ,. ■., Xg"^ , since 
and pg are the largest and the smallest responses. In 

general, the smaller i, the closer to failure the sample Xg 
is. This is shown schematically in Fig. 

Let p S (0,1) be any number such that np is integer. 
By analogy with 0. define the first intermediate failure 
domain Fi as follows: 


Fi={x : g{x) > yl), 


where 


* 


Vi = 



2 


( 20 ) 


( 21 ) 


In other words, Fi is the set of inputs that lead to the 
exceedance of the relaxed threshold 2 /* <y* ■ Note that by 
construction, samples Xg^\ ..., Xg"^^ belong to Fi, while 
Xg"'^'''!),..., Xg") do not. As a consequence, the Direct 
Monte Carlo estimate for the probability of Fi which is 
based on samples Xg^),..., Xg") is automatically equal to 
P, 


1 

nFi)^-Y.lFAx^o^)=P- ( 22 ) 

^ 2=1 


d 

7r(xi,... ,Xd) = (l>[xk), 

fe=i 


( 18 ) 


The value p = 0.1 is often used in the literature, which 
makes Fi a relatively frequent event. Fig. illustrates 
the definition of Fi. 
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FIG. 1: Monte Carlo samples Xq^\ ..., and the failure 
domain F. and are, respectively, the closest to failure 
and the safest samples among Xg^\ ..., Xq"\ 



FIG. 2: The first intermediate failure domain Fi. In this 
schematic illustration, n = 10, p = 0.2, so that there are 
exactly np = 2 Monte Carlo samples in Fi, Xq^\xq^^ G Fi. 


The first intermediate failure domain Fi can be viewed 
as a (very rough) conservative approximation to the tar¬ 
get failure domain F. Since F C Fi, the failure proba¬ 
bility pf can be written as a product: 


PF = P(Fi)P(i^|i^l), 


(23) 


where P(F|Fi) is the conditional probability of F given 
Fi. Therefore, in view of (22), the problem of estimating 
Pf is reduced to estimating the conditional probability 
P(F|Fi). 

In the next stage, instead of generating samples in the 
whole input space (like in DMC), the SS algorithm aims 
to populate Fi. Specifically, the goal is to generate sam¬ 
ples \ ..., from the conditional distribution 


7r(a:|J^i) = 


■K{x)IF^{x) _ / fi {x) 
P(Fi) “ P(Fi) 


W_(t){xk). (24) 




First of all, note that samples Xg^\ ..., Xq"^^ not only be¬ 
long to Fi, but are also distributed according to f(-|Fi). 
To generate the remaining (n — np) samples from f(-|Fi), 
which, in general, is not a trivial task. Subset Simu¬ 
lation uses the so-called Modified Metropolis algorithm 
(MMA). MMA belongs to the class of Markov chain 
Monte Carlo (MCMC ) algorithms (TU [5^, which are 
techniques for sampling from complex probability dis¬ 
tributions that cannot be sampled directly, at least not 
efficiently. MMA is based on the original Metropolis al¬ 
gorithm |23j and specifically tailored for sampling from 
the conditional distributions of the form (24). 


generates another sample x from 7r(-|Fi) as follows: 

1. Generate a “candidate” sample 
For each coordinate k = 1,..., d, 

(a) Sample pk qkiAxk), where qki-\xk), called 
the proposal distribution, is a univariate PDF 
for r]k centered at Xk with the symmetry prop¬ 
erty qkiVklxk) = qkixklVk)- For example, the 
proposal distribution can be a Gaussian PDF 
with mean Xk and variance cr^. 


or it can be a uniform distribution over [xk — 
a ,Xk + a], for some a > 0. 

(b) Gompute the acceptance ratio 


fiixk)' 


(26) 


(c) Define the coordinate of the candidate 
sample by accepting or rejecting pk, 




Pk, with probability min{l,rfe}, 

Xk, with probability 1 — min{l,rfc}. 


(27) 


2. Accept or reject the candidate sample ^ by setting 


X = 


X, if C ^ Fi. 


(28) 


A. Modified Metropolis algorithm 

Let X ~ 7r(-|Fi) be a sample from the conditional dis¬ 
tribution 7r(-|Fi). The Modihed Metropolis algorithm 


The Modified Metropolis algorithm is schematically il¬ 
lustrated in Fig. 1^ 

It can be shown that the sample x generated by MMA 
is indeed distributed according to 7r(-|Fi). If the candi¬ 
date sample ^ is rejected in (28l, then x = x ^ 7r(-|Fi) 
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qk{-\xk) 




FIG. 3: Modified Metropolis algorithm 


and there is nothing to prove. Suppose now that ^ 
is accepted, i = ^, so that the move from x to i is 
a proper transition between two distinct points in Fi. 
Let /(•) denote the PDF of x (the goal is to show that 
/(i) = 7r(i|Fi)). Then 

f{x) = f n{x\Fi)t{x\x)dx, (29) 

JFi 

where t{x\x) is the transition PDF from x to x x. 
According to the first step of MMA, coordinates of i = ^ 
are generated independently, and therefore t{x\x) can be 
expressed as a product. 


t{x\x) = tkiXklXk), 


(30) 


k=l 


where tk{xk\xk) is the transition PDF for the fc*** coordi¬ 
nate Xfc. Combining equations (24), (29), and (30) gives 


fix) = f 

JF 


-fp jx) 

F, nFi) 


n 'Pixk) tk{xk\xk)da 




k^l 


(31) 


P(Ai) 


(j){xk)tk{xk\xk)dx. 


-1 k=i 


The key to the proof of f{x) = 7r(i|Fi) is to demonstrate 
that 4>(xk) and tk(xk\xk) satisfy the so-called detailed 
balance equation, 


(l>{Xk)tkiXk\Xk) = (l){Xk)tkiXk\Xk). 


(32) 


If Xk = Xk, then (32) is trivial. Suppose that Xk ^ Xk, 
that is Xk = Cfc = Fjk in (27). The actual transition 
PDF tk{xk\xk) from Xk to Xk ^ Xk differs from the pro¬ 
posal PDF qk{xk\xk) because the acceptance-rejection 


step (27) is involved. To actually make the move from Xk 
to Xk, one needs not only to generate Xk qk{-\xk), but 
also to accept it with probability minjl, }. There¬ 
fore, 


tk {Xk \xk ) 


qk{xk\xk)min 


0(Xfc) ] 

X^fc)i ’ 


Xk ^ Xk- (33) 


Using ( [^ , the symmetry property of the proposal PDF, 
qk(,Xk\^ = qk{xk\xk), and the identity aminjl,^} = 
&min{l, |} for any a,5 > 0, 


(j){xk)tk{xk\xk) =qkixk\xk)4>{xk) min 

=qkixk\xk)(l){xk) min 
= 4>ixk)tkiXk\Xk), 


(l){Xk) \ 

’Hxk)j 

(l){Xk) \ 


(34) 


and the detailed balance (32) is thus established. The 
rest is a straightforward calculation: 


fix) 


1 

P(Ai) 


d 

(l)ixk)tkixk\xk)dx 

k=l 


1 

P(Ai) 


d 


t{x\x)dx = 7r(x|Fi), 


(35) 


since the transition PDF t{x\x) integrates to 1, and 

IfAx) = 1- 

Remark 2. A mathematically more rigorous proof of the 
Modified Metropolis algorithm is given in 


Remark 3. It is worth mentioning that although the in¬ 
dependence of input variables is crucial for the applica¬ 
bility of MMA, and thus for Subset Simulation, they need 
not be identically distributed. In other words, instead of 
(18), the joint PDF 7r(-) can have a more general form, 
'xix) = Y[k=i^kixk), where nk{-) is the marginal distri¬ 
butions of Xk which is not necessarily Gaussian. In this 
case, the expression for the acceptance ratio in (26) must 
be replaced by rk = 


B. Subset Simulation at higher conditional levels 

Given ~ 7r(-|Fi), it is clear now how to 

generate the remaining (n — np) samples from 7r(-|F’i). 

(i) 

Namely, starting from each Xq , i = 1,..., np, the SS al¬ 
gorithm generates a sequence of (1— ^) new MGMG sam¬ 
ples Xq ^ = Xq I !->■ Xq^{ !->■... !->• xj^^ 1 using the Mod- 
ified Metropolis transition rule described above. Note 
that when Xq) is generated, the previous sample 
is used as an input for the transition rule. The sequence 

(i) (i) (i) 

Xq g, Xq j,..., Xq 1 Is culled u Murkov chain with the 

stationary distribution 7r(-|Fi), and XqI = Xq^ is often 
referred to as the “seed” of the Markov chain. 
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FIG. 4: MCMC samples generated by the Modified Metropo¬ 
lis algorithm at the first conditional level of Subset Simula¬ 
tion. 


FIG. 5: The second intermediate failure domain F 2 . In this 
schematic illustration, n = 10, p = 0.2, so that there are exactly 
np = 2 MCMC samples in F 2 , E F 2 . 


To simplify the notation, denote samples 
i-i • ■ ■! The subscript 1 in¬ 


dicates that the MCMC samples ~ 7 r(-|Fi) 

are generated at the first conditional level of the SS 
algorithm. These conditional samples are schematically 
shown in Fig. Also assume that the corresponding 

j. (1) ! (1)\ O' 

system responses y\ = g{x\ ,yl 


= g{x^i^) are 


arranged in the decreasing order, i.e. y^'^ > .. • > ■ 

If the failure event F is rare enough, that is if pf is 
sufficiently small, then it is very likely that none of the 
samples belongs to F, i.e. y^^ < y* for all 

i = l,...,n. Nevertheless, these MCMC samples can 
be used in the similar way the Monte Carlo samples 


.(1) 


(n) 1 

Xn were used. 


By analogy with (20), define the second intermediate 


failure domain F 2 as follows: 


F 2 = {x : g{x) > 2 / 2 }, 


where 


2/2 = 


Vi 


{np) 


■Vi 


(np+1) 


(36) 


(37) 


Note that y^ > y* since y^^ > y^ for alH = 1,... ,n. 
This means that F C F 2 C Fi, and therefore, F 2 can be 
viewed as a conservative approximation to F which is still 
rough, yet more accurate than Fi. Fig. [^illustrates the 
definition of F 2 . By construction, samples x^^\ ..., 
belong to F 2 , while ..., x^^ do not. As a result, 

the estimate for the conditional probability of F 2 given 
Fi which is based on samples ,... ~ 7 r(-|Fi) is 

automatically equal to p. 


1 "■ 

P(F2|Fi)«-5]/f,(4^)=P. 


(38) 


Since F C F 2 C Fi, the conditional probability 


P(F|Fi) that appears in (23) can be expressed as a prod¬ 
uct: 


P(F|Fi) =P(F2 |Fi)P(F|F2). 


(39) 


Combining (23) and (39) gives the following expression 


for the failure probability: 


Pf=P(Fi)P(F 2 |Fi)P(F|F 2 ). 


(40) 


Thus, in view of (22) and (38 1 , the problem of estimating 


Pf is now reduced to estimating the conditional proba¬ 
bility P(F|F 2 ). 

In the next step, as one may have already guessed, the 
Subset Simulation algorithm: populates F 2 by generat¬ 
ing MCMC samples x^\ ..., from 7 r(-|F 2 ) using the 
Modified Metropolis algorithm; defines the third inter¬ 
mediate failure domain F 3 C F 2 such that P(F 3 |F 2 ) « 
^ /F 3 (a; 2 *^) = p ; and reduces the original prob¬ 

lem of estimating the failure probability pf to estimat¬ 
ing the conditional probability P(F|F 3 ) by represent¬ 
ing Pf = P(Fi)P(F 2 |Fi)P(F 3 |F 2 )P(F|F 3 ). The algo¬ 
rithm proceeds in this way until the target failure do¬ 
main F has been sufficiently sampled so that the condi¬ 
tional probability P(F|Ff) can be accurately estimated 
by i X:r=i where Fp is the intermediate fail¬ 


ure domain, and x\ 


( 1 ) 




7 r(-|FF) are the MCMC 


samples generated at the conditional level. Subset 
Simulation can thus be viewed as a method that decom¬ 
poses the rare failure event F into a sequence of progres¬ 
sively “less-rare” nested events, F C Ff C ... C Fi, 
where all intermediate failure events Fi,..., Fp are con¬ 
structed adaptively by appropriately relaxing the value 
of the critical threshold yl < ... < y^ < y*. 











C. stopping criterion 


In what follows, the stopping criterion for Subset Simu¬ 
lation is described in detail. Let 71^(0 denote the number 
of failure samples at the level, that is 

n 

^f{1) = (41) 

2=1 


where ^ 7r(-|J^;). Since is a rare event, 

it is very likely that np(l) = 0 for the first few condi¬ 
tional levels. As I gets larger, however, np{l) starts in¬ 
creasing since Fi , which approximates F “from above”, 
shrinks closer to F. In general, np{l) > np{l — 1), since 
F C F, C Fi_i and the np closest to F samples among 
, ■ ■ ■, xl_\ are present among x\ \ ... ,xl . At con¬ 
ditional level I, the failure probability pf is expressed as 
a product. 


pp = P(Fi)P(F2|Fi) . ..F{Fi\Fi_i)F{F\Fi). (42) 


Furthermore, the adaptive choice of intermediate critical 
thresholds pi,... ,y* guarantees that the first I factors in 
(42) approximately equal to p, and, thus. 


The described stopping criterion guarantees that the 
estimated values of all factors in the factorization pp = 
P(Fi)P(A 2 |F’i) ...P(Af|Ff_i)P(F"|F"f) are not smaller 
than p. If p is relatively large {p = 0.1 is often used 
in applications), then it is likely that the estimates 
P(Fi) « p, P(A 2 |Fi) « p,...,F{Fl\Fl_i) « p, and 
P(F|Ff) « (^ p) are accurate even when the sam¬ 

ple size n is relatively small. As a result, the SS estimate 
(45) is also accurate in this case. This provides an intu¬ 
itive explanation as to why Subset Simulation is efficient 
in estimating small probabilities of rare events. For a de¬ 
tailed discussion of error estimation for the SS method 
the reader is referred to [3]. 


D. Implementation details 

In the rest of this section, the implementation details 
of Subset Simulation are discussed. The SS algorithm 
has two essential components that affect its efficiency: 
the parameter p and the set of univariate proposal PDFs 
{qk}, k = l,...,d. 


1. Level probability 


pp^p^-FiFlFi). 


(43) 


Since there are exactly np{l) failure samples at the 
leve l, the estimate of the last conditional probability in 
(42) which is based on samples ..., ^ 7r('|^i) is 

given by 


F{F\Fi) 


n 

-Y.Fi 

n A—^ 


) = 


np{l) 


(44) 


i=l 


If np{l) is sufficiently large, i.e. the conditional event 


(F"|Fi) is not rare, then the estimate (44) is fairly accu¬ 


rate. This leads to the following stopping criterion: 

• If > p, i.e. there are at least np failure sam¬ 
ples among ..., then Subset Simulation 
stops: the current conditional level I becomes the 
last level, L = I, and the failure probability esti¬ 
mate derived from (43) and (|44| is 


PF « Pf = 


(45) 


The parameter p, called the level probability in [5] and 
the conditional failure probability in |29j . governs how 
many intermediate failure domains F] are needed to reach 
the target failure domain F. As it follows form (45), a 


small value of p leads to a fewer total number of con¬ 
ditional levels L. But at the same time, it results in a 
large number of samples n needed at each conditional 
level I for accurate determination of F] (i.e. determina¬ 
tion of y*) that satisfies ^ X)r=i(^i-\) = P- the 
extreme case when p < pf, no levels are needed, F = 0, 
and Subset Simulation reduces to the Direct Monte Carlo 
method. On the other hand, increasing the value of p will 
mean that fewer samples are needed at each conditional 
level, but it will increase the total number of levels L. 
The choice of the level probability p is thus a tradeoff 
between the total number of level L and the number of 
samples n at each level. In the original paper [T], it has 
been found that the value p = O.I yields good efficiency. 
The latter studies [3 [29], where the c.o.v. of the SS es¬ 
timate jfp has been analyzed, confirmed that p = O.I is 
a nearly optimal value of the level probability. 


• If < p, i.e. there are less than np fail¬ 

ure samples among ..., then the algo¬ 
rithm proceeds by defining the next intermediate 
failure domain = {a; : g(x) > y*_f_i}, where 

Pi+i = ( 2 //”^^ + and expressing P(F|F;) 

as a product P(F|F)) = P(F;+i|F;)P(F|Fi+i) « 
p-P(F|Ti+i). 


2. Proposal distributions 

The efficiency and accuracy of Subset Simulation also 
depends on the set of univariate proposal PDFs {qt}, 
k = 1,... ,d that are used within the Modified Metropo¬ 
lis algorithm for sampling from the conditional distri¬ 
butions 7r(-|F;). To see this, note that in contract to 
the Monte Carlo samples a;g^\ ..., x^f^^ ~ 7r(-) which are 
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i.i.d., the MCMC samples xl 


( 1 ) 




7r(-|Fi) are 


not independent for / > 1, since the MMA transition 


rule uses ~ to generate ~ Tr{-\Fi). 

This means that although these MCMC samples can 
be used for statistical averaging as if they were i.i.d., 
the efficiency of the averaging is reduced if compared 
with the i.i.d. case [H]. Namely, the more correlated 
..., x\^'^ are, the slower is the convergence of the es¬ 
timate P{Fi+i\Fi) « and, therefore, 

the less efficient it is. The correlation between samples 
..., x\^'^ is due to proposal PDFs {qk}, which govern 
the generation of the next sample from the current 

one Hence, the choice of {qk} is very important. 

It was observed in [T] that the efficiency of MMA is 
not sensitive to the type of the proposal PDFs (Gaus¬ 
sian, uniform, etc), however, it strongly depends on their 
spread (variance). Both small and large spreads tend 
to increase the correlation between successive samples. 


d*+i) 


Large spreads may reduce the acceptance rate in (281, in¬ 
creasing the number of repeated MCMC samples. Small 
spreads, on the contrary, may lead to a reasonably high 
acceptance rate, but still produce very correlated sam¬ 
ples due to their close proximity. As a rule of thumb, 
the spread of qk, k = 1,... ,d, can be taken of the same 
order as the spread of the corresponding marginal PDF 
'^k 0. For example, if tt is given by ( [l^ , so that all 


marginal PDFs are standard Gaussian, Trk{x) = (j){x), 
then all proposal PDFs can also be Gaussian with unit 
variance, qk{x\xk) = ^(x — Xk). This choice is found to 
give a balance between efficiency and robustness. 

The spread of proposal PDFs can also be chosen adap¬ 
tively. In |29j . where the problem of optimal scaling for 
the Modified Metropolis algorithm was studied in more 
detail, the following nearly optimal scaling strategy was 
proposed: at each conditional level, select the spread 
such that the the corresponding acceptance rate in (28) 


is between 30% and 50%. In general, finding the optimal 
spread of proposal distributions is problem specific and 
a highly non-trivial task not only for MMA, but also for 
almost all MGMG algorithms. 


VI. ILLUSTRATIVE EXAMPLES 


To illustrate Subset Simulation and to demonstrate its 
efficiency in estimating small probabilities of rare failure 
events, two examples are considered in this section. As it 
has been discussed in Section |llj in reliability problems, 
the dimension d of the input space X is usually very 
large. In spite of this, for visualization and educational 
purposes, a linear reliability problem in two dimensions 
(d = 2) is first considered in Section VIA A more realis¬ 
tic high-dimensional example (d = 10'^) is considered in 
the subsequent Section VI B[ 


A. Subset Simulation in 2-D 

Suppose that d = 2, i.e. the response variable y de¬ 
pends only on two input variables xi and X 2 . Consider a 
linear performance function. 


g{Xi,X2) =Xi+ X2, 


(46) 


where Xi and X 2 are independent standard Gaussian, 
Xi ~ i = 1,2. The failure domain F is then 

a half-plane defined by 


F = {{xi,X 2) : xi+X2>y*}. 


(47) 


In this example, the failure probability pp can be cal¬ 
culated analytically. Indeed, since xi + X 2 ^ A/’(0,2), 
and, therefore, ^ Af{0, 1), 


Pf = P(a:i +X 2 > y*) = ¥ [ ^ 


y 




(48) 


where $ is the standard Gaussian CDF. This expression 
for the failure probability can be used as a check on the 
SS estimate. Moreover, expressing y* in terms of pp. 


y 


= V2^-\1-Pf), 


(49) 


allows to solve the inverse problem, namely, to formu¬ 
late a linear reliability problem with a given value of the 
failure probability. Suppose that pp = 10“^° is the tar¬ 
get value. Then the corresponding value of the critical 
threshold is y* « 9. 

Subset Simulation were used to estimate the failure 


probability of the rare event (47) with y* = 9. The pa¬ 
rameters of the algorithm were chosen as follows: the 
level probability p = 0.1, the proposal PDFs qk{x\xk) = 
4>{x — Xk), and the sample size n = 10^ per each level. 
This implementation of SS led to L = 9 conditional 
levels, making the total number of generated samples 
N = n + L{n — np) = 9.1 x 10^. The obtained SS es¬ 
timate is Pp = 1.58 X 10“^° which is quite close to the 
true value pp = 10“^°. Note that, in this example, it 
is hopeless to obtain an accurate estimate by the Direct 


Monte Carlo method since the DMC estimate (12) based 
on TV = 9.1 X 10^ samples is effectively zero: the rare 
event F is too rare. 

Fig. I^shows the samples generated by the SS method. 
The dashed lines represent the boundaries of intermedi¬ 
ate failure domains Fi, I = 1,... ,L = 9. The solid line 
is the boundary of the target failure domain F. This 
illustrates how Subset Simulation pushes Monte Carlo 
samples (red) towards the failure region. 


B. Subset Simulation in High Dimensions 

It is straightforward to generalize the low-dimensional 
example considered in the previous section to high di- 
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FIG. 6: Samples generated by Subset Simulation: red samples 
are Monte Carlo samples generated at the 0**“ unconditional 
level, purple sample are MCMC sample generated at the I'’* 
conditional level, etc. The dashed lines represent the bound¬ 
aries of intermediate failure domains Fi, I = 1,...,L = 9. 
The solid line is the boundary of the target failure domain F. 
[Example 6.1]. 


mensions. Consider a linear performance function 

d 

9{x)='^Xi, (50) 

i=l 

where Xi,... ,Xd are i.i.d. standard Gaussian. The fail¬ 
ure domain is then a half-space defined by 


d 

F = {x : '^Xi> y*}. (51) 

i=l 


In this example, d = 10^ is considered, hence the input 
space is indeed high-dimensional. As before, the 

failure probability can be calculated analytically: 


Pf 



>y 



\ Vd Vd) ^^-2) 


This expression will be used as a check on the SS esti¬ 
mate. 

First, consider the following range of values for the 
critical threshold, y* G [0,200]. Fig. plots y* ver¬ 
sus Pf- The solid red curve corresponds to the sample 
mean of the SS estimates p^p which is based on 100 in¬ 
dependent runs of Subset Simulation. The two dashed 
red curves correspond to the sample mean ± one sample 
standard deviation. The SS parameters were set as fol¬ 
lows: the level probability p = 0.1, the proposal PDFs 


FIG. 7: Critical threshold y* versus the failure probability 
Pf- [Example 6.2]. 


qk{x\xk) = (j){x — Xk), and the sample size n = 3x 10^ per 
each level. The solid blue curve (which almost coincides 
with the solid red curve) corresponds to the true values 
oi Pf computed from (52). The dark green curves cor¬ 
respond to Direct Monte Carlo: the solid curve is the 
sample mean (based on 100 independent runs) of the 
DMC estimates (12), and the two dashed curves 

are the sample mean ± one sample standard deviation. 
The total number of samples N used in DMC equals to 
the average (based on 100 runs) total number of sam¬ 
ples used in SS. Finally, the dashed light green curves 
show the theoretical performance of Direct Monte Carlo, 
namely, they correspond to the true value of pf (52) ± 
one theoretical standard deviation obtained from (14). 
The bottom panel of Fig. shows the zoomed in region 
that corresponds to the values y* G [100,160] of the crit¬ 
ical threshold. Note that for relatively large values of the 
failure probability, pf < 10“^, both DMC and SS pro¬ 
duce accurate estimates of pf- For smaller values how¬ 
ever, Pf < 10“®, the DMC estimate starts to degenerate, 
while SS still accurately estimates pf- This can be seen 
especially well in the bottom panel of the figure. 

The performances of Subset Simulation and Direct 
Monte Carlo can be also compared in terms of the co¬ 
efficient of variation of the estimates p^ and . This 
comparison is shown in Fig. The red and dark green 
curves represent the sample c.o.v. for SS and DMC, re¬ 
spectively. The light green curve is the theoretical c.o.v. 
of p™'^ given by (15). When the critical threshold is rel¬ 
atively small y* < 60, the performances of SS and DMC 
are comparable. As y* gets large, the c.o.v. of 
starts to grow much faster than that of pp*. In other 
words, SS starts to outperform DMC, and the larger y*, 
i.e. the more rare the failure event, the more significant 
the outperformance is. 
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FIG. 8: Critical threshold y* versus the c.o.v. [Example 6.2]. 



Failure probability 

FIG. 9: Failure probability versus the total number of sam¬ 
ples. [Example 6.2]. 



Independent runs 


FIG. 10: Performance of Subset Simulation for 100 independent 
runs. The critical threshold is y* = 200, the corresponding true 
value of the failure probability is pp = 1.27 X lO”^'^. [Example 
6 . 2 ]. 



FIG. 11: System responses > ... > n = 3 X 10®, for all 
levels, Z = 0,..., L = 9, for a fixed simulation run. [Example 6.2]. 


The average total number of samples used in Sub¬ 
set Simulation versus the corresponding values of failure 
probability is shown in the top panel of Fig.|^ The stair¬ 
case nature of the plot is due to the fact that every time 
when Pp crosses the value by decreasing from p^ + e 
to p^ — e, an additional conditional level is required. In 
this example, p = 0.1 is used, that why the jumps occur 
at Pp = 10“^, k = 1,2,_ The jumps are more pro¬ 

nounced for larger values of pp, where the SS estimate 
is more accurate. For smaller values ofpp, where the SS 
estimate is less accurate, the jumps are more smoothed 
out by averaging over independent runs. 

In Fig. in where the c.o.v’s of SS and DMC are com¬ 
pared, the total numbers of samples (computational ef¬ 
forts) used in the two methods are the same. The natural 
question is then the following: by how much should the 
total number of samples N used in DMC be increased 
to achieve the same c.o.v as in SS (so that the green 
curve in Fig. m coincides with the red curve)? The an¬ 
swer is given in the bottom panel of Fig. For example, 
if Pp = 10“^°, then N = 10^°, while the computational 
effort of SS is less than 10^ samples. 

Simulation results presented in Figures |7|8| and 


clearly indicate that (a) Subset Simulation produces a 
relatively accurate estimate of the failure probability, 
and (b) Subset Simulation drastically outperforms Di¬ 
rect Monte Carlo when estimating probabilities of rare 
events. 


Let us now focus on a specific value of the critical 
threshold, y* = 200, which corresponds to a very rare 
failure event (51) with probability pp = 1.27 x 10“^°. 
Fig. demonstrates the performance of Subset Simu¬ 
lation for 100 independent runs. The top panel shows 
the obtained SS estimate p^p for each run. Although 
Pp varies signihcantly (its c.o.v. is S{p^) = 0.74), its 
mean value pp = 1.18 x 10“^° (dashed red line) is close 
to the true value of the failure probability (dashed blue 
line). The bottom panel shows the total number of sam¬ 
ples used in SS in each run. It is needless to say that 
the DMC estimate based on ~ 3 x 10"*^ samples would 
almost certainly be zero. 


Fig. shows the system responses > .. • > 
n = 3 X10®, for all levels, Z = 0,..., L = 9, for a fixed sim¬ 
ulation run. As expected, for the first few levels (6 levels 
in this case), the number of failure samples np{l), i.e. 
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Conditional level 

FIG. 12: Intermediate critical thresholds L = 9, at 

different conditional levels in a fixed simulation run. [Example 
6 . 2 ]. 


samples x'f'^ with = g{xf'’) > y*, is zero. As Sub¬ 


set Simulation starts pushing the samples towards the 
failure domain, np^l) starts increasing with npid) = 3, 
^f(7) = 6, npiS) = 59, and, finally, npid) = 582, af¬ 
ter which the algorithm stopped since nF{9)/n = 0.194 
which is large than p = 0.1. Finally, Fig. [T^ plots the in¬ 
termediate (relaxed) critical thresholds j/)",..., at dif¬ 
ferent levels obtained in a fixed simulation run. 


VII. MATLAB CODE 

This section contains the MATLAB code for the exam¬ 
ples considered in Section [Vl| For educational purposes, 
the code was written as readable as possible with numer¬ 
ous comments. As a result of this approach, the efficiency 
of the code was unavoidably scarified. This code is also 
available online at http: //arxiv. org/, 
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% Subset Simulation for Liner Reliability Problem 
% Performance function: g(x)=xl+ ... +xd 
% Input variables xl ,..., xd are i.i.d. N(0,1) 

% Written by K.M. Zuev, Institute of Risk & Uncertainty, Uni of Liverpool 


clear; 

d=1000; 

YF=200; 

pF=l—normcdf {YF/sqrt (d) ) ; 
n=3000; 

P=0.1; 
nc=n*p; 
ns=(1—p)/p; 


% dimension of the input space 
% critical threshold (failure <=> g(x)>YF) 
% true value of the failure probability 
% number of samples per level 
% level probability 
% number of Markov chains 
% number of states in each chain 


L=0; 

x=randn(d,n); 

nF=0; 

for i=l:n 

y (i)=sum{x{: , i) ) ; 
if y (i)>YF 
nF=nF+l; 

end 


% current (unconditional) level 
% Monte Carlo samples 
% number of failure samples 

% system response y=g(x) 

% y(i)>YF <=> x(:,i) is a failure sample 


end 

while nF(L+l)/n<p % stopping criterion 

L=L+1; % next conditional lelvel is needed 

[y(L,:),ind]=sort(y(L, descend ') ; % renumbered responses 
X(:,:,L)=x(:,ind(:),L); % renumbered samples 

Y (L) = (y (L,nc) +y (L,nc+1) )/2; % L^th intermediate threshold 

z (:, :, 1) =x (:, 1: nc, L) ; % Markov chain ’’ seeds " 


% Modified Metropolis algorithm for sampling from pi(x|F_L) 
for j=l:nc 

for m=l:ns 

% Step 1 : 
for k=l:d 

a=z(k,j,m)+randn; % Step 1(a) 

r=min(1,normpdf(a)/normpdf(z(k,j,m))); % Step 1(b) 

% Step 1 (c) : 
if rand<r 
q(k)=a; 

else 

q(k)=z(k,j,m); 

end 

end 

% Step 2 : 


45 
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46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 


if sum(q)>Y(L) % q belongs to F_L 
z , j,m+l) =q; 

else 

z(:,j,m+1)= z(:,j,m); 

end 

end 

end 

for j=l:no 

for m=l:ns+1 

x(:, (j—l)*(ns+l)+m,L+l)=z(:,j,m); % samples from pi(x|F_L) 

end 

end 

clear z; 


nF(L+1)=0; 
for i=l:n 

y(L+1, i)=sum(x(:,i,L+1)); % system response y=g(x) 
if y(L+l,i)>YF % then x(:,i,L+l) is a failure sample 

nF(L+1)=nF(L+1)+1; % number of failure samples at level L+1 

end 

end 

end 

pF_SS=p^ (L)*nF(L+1)/n; % SS estimate 

N=n+n*(1—p)*(L); % total number of samples 


VIII. SUMMARY 

In this paper, a detailed exposition of Subset Simula¬ 
tion, an advanced stochastic simulation method for esti¬ 
mation of small probabilities of rare events, is provided 
at introductory level. A simple step-by-step derivation 
of Subset Simulation is given, and important implemen¬ 
tation details are discussed. The method is illustrated 
with a few intuitive examples. 

After the original paper [1] was published, various 
modifications of SS were proposed: SS with splitting [^, 
Hybrid SS [B] , and Two-Stage SS [T3] , to name but a few. 
It is important to highlight, however, that none of these 
modifications offer a drastic improvement over the orig¬ 
inal algorithm. A Bayesian analog of SS was developed 
in |29| . For further reading on Subset Simulation and its 
applications, a fundamental and very accessible mono¬ 
graph [3] is strongly recommended, where the method 


is presented from the CCDF (complementary cumulative 
distribution function) perspective and where the error 
estimation is discussed in detail. 

Also, it is important to emphasize that Subset Sim¬ 
ulation provides an efficient solution for general relia¬ 
bility problems without using any specific information 
about the dynamic system other than an inputoutput 
model. This independence of a systems inherent proper¬ 
ties makes Subset Simulation potentially useful for appli¬ 
cations in different areas of science and engineering. 

As a final remark, it is a pleasure to thank Professor 
Siu-Kui Au whose comments on the first draft of the pa¬ 
per were very helpful. Professor James Beck, who gener¬ 
ously shared his knowledge of and experience with Subset 
Simulation, and Professor Francis Bonahon for general 
support and for creating a nice atmosphere at the De¬ 
partment of Mathematics of the University of Southern 
California, where the author started this work. 
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